//  Denis Yurin (c) 2012

#ifndef RADIAL_PROFILE_FINDER_H
#define RADIAL_PROFILE_FINDER_H

#include <iostream>
#include <string.h> // memset
#include <cstdlib> // std::system("...");

#include "Vector/vector.03.h"
#include "GadgetParticleType/gpt.01.h"
#include "utils/utils.02.h"


class RPF {
	
	public:
		
		RPF():
		r(0),
		n_bin(0),
		nc_bin(0),
		m_bin(0),
		mc_bin(0),
		vt_bin(0),
 		rho_bin(0),
 		rhoc_bin(0),
		G(43.0071),
		H0(100),
		a(1),
		r1(0),
		rho_crit(3*H0*H0/(8*pi*G)),
		Rvir(1),
		logflag(false)
		{}
		
		double *r;	// radius
		
		uint *n_bin;	// shell number bins
		uint *nc_bin;	// comulative number bins
		double *m_bin; 	// shell mass bins
		double *mc_bin; 	// comulative mass bins
		
		double *vt_bin;	// tangential velocity bins: vt_bin[0] at 0-r0, vt_bin[bin_n+1] at r2-infinity
		
		
 		double *rho_bin;	// shell density bins
		double *rhoc_bin; // enclosed density bins
		
		double G;
		double H0;
		double a;
		double rho_crit;
		double Rvir;
		double  r1,  r2;
		double _r1, _r2;
		uint bin_n; 
		double bin_size;
		static const double pi = 3.141592654;
		bool logflag;
		
		void find(uint N, uint ofs, P *p, Vector<double> C, double r1, double r2, uint bin_n, double a=1){
			
			float *m = new float[N];
			Vector<float> *R = new Vector<float>[N];
			Vector<float> *V = new Vector<float>[N];
			
			for (uint i=0; i<N; i++) {
				m[i] = p[ofs+i].Mass;
				R[i] = *(Vector<float>*)p[ofs+i].Pos;
				V[i] = *(Vector<float>*)p[ofs+i].Vel;
			}
			
			find(N, m, R, V, C, r1,r2, bin_n, a);
			
			delete [] R, V;
			delete [] m;
		}
		
		template <typename T>		
		void find(uint N, T *m, Vector<T> *R, Vector<T> *V, Vector<double> C, double r1, double r2, uint bin_n, double a=1){
			
			// please note that bins are calculated in proper coordinates (r1, r2) 
			
			this->bin_n = bin_n;
			this->a = a;
			this->r1 = r1;
 			this->r2 = r2;
			this->_r1 = r1/a;
 			this->_r2 = r2/a;
	

			bin_size = logflag ? log10(r2/r1) / bin_n : (r2 - r1) / bin_n;
			
			// bin's coordinates
			if (r) delete [] r;
			r = new double[bin_n+1]; 
			for (uint i=0; i<=bin_n; i++)
			r[i] = logflag ? pow(10, log10(r1) + i*bin_size) : r1 + i*bin_size;
			
			// Shell Mass & Number BINs
			if (m_bin) delete [] m_bin;
			m_bin = new double[bin_n+2];
			memset(m_bin, 0, bin_n*sizeof(double));
			
			if (n_bin) delete [] n_bin;
			n_bin = new uint[bin_n+2];
			memset(n_bin, 0, bin_n*sizeof(uint));
			
			if (V) {
				if (vt_bin) delete [] vt_bin;
				vt_bin = new double[bin_n+2]; 
				memset(vt_bin, 0, (bin_n+2)*sizeof(double));
			}
			
			for (uint i=0; i<N; i++) {
				
				Vector<T> DR = R[i]-C;
				double dr= a*DR.module();  // NB: Here we acctually do converttion from comove to proper
				
				uint bin_id = logflag ? log10(dr/r1) / bin_size : (dr-r1) / bin_size;
			
				// Logic: [0] E 0-r1, [1-N] E r1-r2  [N+1] E r2-inf
				if (bin_id==-1) bin_id=0; else bin_id++;
				if (bin_id>bin_n) bin_id = bin_n+1;
				
				m_bin[bin_id]+= m[i];
				n_bin[bin_id]++;
				
 				if (V) 
				vt_bin[bin_id]+=pow(V[i].module(),2);
				//if (V) rcr_bin[bin_id]+=(V[i]-(V[i]*DR.ort())*DR.ort()).module();
				
				
			}
			
			if (V)
			for (uint i=0; i<=bin_n+1; i++)
			vt_bin[i]= n_bin[i] ? sqrt(vt_bin[i]/n_bin[i]) : -1;
			
			
			
			
			// Comulative Mass & Number bins
			if (mc_bin) delete [] mc_bin;
			mc_bin = new double[bin_n+2]; 
			
			if (nc_bin) delete [] nc_bin; 
			nc_bin = new uint[bin_n+2]; 
			
			mc_bin[0]=m_bin[0];
			nc_bin[0]=n_bin[0];
			for (uint i=1; i<=bin_n+1; i++) {
				mc_bin[i]=mc_bin[i-1]+m_bin[i];
				nc_bin[i]=nc_bin[i-1]+n_bin[i];
			}
			
			// Spherical Density bins
			if (rhoc_bin) delete [] rhoc_bin;
			rhoc_bin = new double[bin_n+1];
			for (uint i=0; i<=bin_n; i++) {
				double  Vi = (4*pi/3)*pow(r[i],3);
				rhoc_bin[i] = Vi ? mc_bin[i] / Vi : -1;
			}
			
			// Shell Density bins
			if (rho_bin) delete [] rho_bin;
			rho_bin = new double[bin_n+1];
			*rho_bin = *rhoc_bin;
			for (uint i=1; i<=bin_n; i++)
			rho_bin[i] =  m_bin[i] / ( 4.0*pi/3.0 * ( pow(r[i],3.0) - pow(r[i-1],3.0) ) );

		}


		void print(){
			
			printf("\nRadialProfileFinder::print():\n\n");
			
			logflag ? printf("logflag = true\n") : printf("logflag = false\n");
			printf("rho_crit = %g\n", rho_crit);
			

			for (uint i=1; i<=bin_n; i++)
				printf("%2u % 7.3e    \t %7u  % 7.3e % 7.3e % 7.3e % 7.3e % 8.3f\n", 
						i, r[i], 
						n_bin[i], m_bin[i], mc_bin[i], rho_bin[i], rhoc_bin[i], vt_bin[i]);
			
		}
		
		void plot(const char* filename){
			/*
			char cstr[10000];			
			
			// 1. Gnuplot filename.rp
			sprintf(cstr, "%s.radialProfiles.dat", filename);
			FILE *file = fopen(cstr,"w");
			for (uint i=1; i<=bin_n; i++) {
				double rmi = logflag ? pow(10,log10(r[i])-0.5*bin_size) : 0.5*(r[i+1]-r[i]);
				fprintf(file, "%10.7e %10.7e %10u %10.7e %10.7e %10.7e %10.7e %10.7e\n", r[i], rmi, n_bin[i], m_bin[i], mc_bin[i], rho_bin[i], rhoc_bin[i], rc_bin[i]);
			
			}
			fclose(file);
			
			// 2. Gnuplot filename.rp.plt 
			std::string plt;
			
			char cstr1[] = {
				"set term postscript landscape enhanced color 'Text' 12\n"
				"set output '%s.radialProfiles.ps'\n"
				"set key top left\n"
				"set logscale xy\n"
				"set grid xtics mxtics ytics mytics lt 1 lc rgb '#CCCCCC' lw 0.5 \n"
				"\n"
			};
			
			sprintf(cstr, cstr1, filename);
			plt+=cstr;
			
			// Density
			char cstr2[] = {

				"set title 'Enclosed Particle Number: N(<r)'\n"
				"plot '%s.radialProfiles.dat' u 1:3	t '' 	w l lw 4 lt 1 lc 1\n"
				
				"set title 'Shell Mass: M(r)'\n"
				"plot '%s.radialProfiles.dat' u 1:4	t '' 	w l lw 4 lt 1 lc 1\n"
				
				"set title 'Enclosed Mass: M(<r)'\n"
				"plot '%s.radialProfiles.dat' u 1:5	t '' 	w l lw 4 lt 1 lc 1\n"
				
				"set title 'Shell Density: rho(r)'\n"
				"plot '%s.radialProfiles.dat' u 2:6	t '' 	w l lw 4 lt 1 lc 1\n"
				
				"set title 'Enclosed Density: rho(<r)'\n"
				"plot '%s.radialProfiles.dat' u 1:7	t '' 	w l lw 4 lt 1 lc 1\n"
				
				"set title 'Rotation Curve: V(r)'\n"
				"plot '%s.radialProfiles.dat' u 1:8	t '' 	w l lw 4 lt 1 lc 1\n"
				
			};
			
			sprintf(cstr, cstr2, filename, filename, filename, filename, filename, filename);
			plt+=cstr;
			
			FILE *file_plt = fopen((string(filename)+string(".radialProfiles.plt")).c_str(),"w");
			fprintf(file_plt, "%s", plt.c_str());
			fclose(file_plt);
			
			// 3. run gnuplot filename.rp.plt 
			sprintf(cstr, "gnuplot %s.radialProfiles.plt", filename);
			std::system(cstr);
			
			// 4. run ocular filename.rp.ps 
			sprintf(cstr, "okular %s.radialProfiles.ps", filename);
			std::system(cstr);
			*/
		}
		

};



#endif